
using System;

namespace PVTLibrary.MathEx
{
    public class Brent
    {

        public funcdelegate _fc;

        public delegate double funcdelegate(double x, object otherargs);


        public void DefineFuncDelegate(funcdelegate fg)
        {
            _fc = fg;
        }


        double func(double x , object otherargs ) 
        {
            return _fc.Invoke(x, otherargs);


        }


        public double BrentOpt(double minval, double maxval, int n, double tol, int itmax, object otherargs)
        {
            double x_inf, x_sup, y_inf, y, delta_x;

            x_inf = minval;
            x_sup = maxval;

            delta_x = (x_sup - x_inf)/n;

            do
            {
                y = func(x_inf, otherargs);
                x_inf = x_inf + delta_x;
                y_inf = func(x_inf, otherargs);
            } while (y*y_inf < 0 || x_inf > x_sup);

            x_sup = x_inf - delta_x;

            double  aaa, bbb, ccc, ddd, eee, min11, min22, faa, fbb, fcc, ppp, qqq, rrr, sss, tol11, xmm;
            int ITMAX2 = itmax;
            int iter2;

            aaa = x_inf;
            bbb = x_sup;
            ccc = x_sup;

            faa = func(x_inf, otherargs);
            fbb = func(x_sup, otherargs);
            fcc = func(x_sup, otherargs);

            iter2 = 0;
            do
            {                
                if( (fbb > 0 && fcc > 0) || (fbb < 0 && fcc < 0) )
                {
                    ccc = aaa;
                    fcc = faa;
                    ddd = bbb - aaa;
                    eee = ddd;
                }

                if(Math.Abs(fcc) < Math.Abs(fbb) )
                {
                    aaa = bbb;
                    bbb = ccc;
                    ccc = aaa;
                    faa = fbb;
                    fbb = fcc;
                    fcc = faa;
                }
                tol11 = tol;
                xmm = 0.5 * (ccc - bbb);

                if( (Math.Abs(xmm) <= tol11) || (fbb == 0)) goto Final3;
                
                if( Math.Abs(fbb) < tol11)  goto Final3;
                
                if ( Math.Abs(eee) >= tol11 && Math.Abs(faa) > Math.Abs(fbb)) 
                {
                    sss = fbb / faa;
                    if (aaa == ccc)
                    {
                        ppp = 2*xmm*sss;
                        qqq = 1 - sss;
                    }
                    else
                    {
                        qqq = faa/fcc;
                        rrr = fbb/fcc;
                        ppp = sss*(2*xmm*qqq*(qqq - rrr) - (bbb - aaa)*(rrr - 1));
                        qqq = (qqq - 1)*(rrr - 1)*(sss - 1);
                    }

                    if (ppp > 0) qqq = -qqq;

                    ppp = Math.Abs(ppp);
                    min11 = 3*xmm*qqq - Math.Abs(tol11*qqq);
                    min22 = Math.Abs(eee*qqq);
                    double tvar2;

                    if (min11 < min22) tvar2 = min11;
                    if (min11 > min22) tvar2 = min22;
                    if( 2 * ppp < tvar2 )
                    {
                        eee = ddd;
                        ddd = ppp/qqq;
                    }
                    else
                    {
                        ddd = xmm;
                        eee = ddd;
                    }
                }
                else
                {
                    ddd = xmm;
                    eee = ddd;
                }
                aaa = bbb;
                faa = fbb;
                if (Math.Abs(ddd) > tol11)
                {
                    bbb += ddd;
                }
                else
                {
                    bbb += Math.Sign(xmm)*tol11;
                }
                fbb = func(bbb, otherargs);
                iter2 += 1;
            }
            while(iter2 == ITMAX2);

            return bbb;

Final3:
            return bbb;

        }
    }
}